Bioactivity of dihydropyrimidinone derivatives as inhibitors of cyclooxygenase-2 (COX-2): an in silico approach

Cyclooxygenase-2 (COX-2) is an enzyme involved in inflammation. The overexpression of COX-2 causes chronic inflammation, which can be prevented by COX-2 inhibitors. Generally, COX-2 inhibitors possess a carboxyl group and an aromatic ring in their molecular structure. These moieties are involved in the interaction with the active site of COX-2, thus playing a pivotal role in the inhibitory activity. Regarding the requisite molecular structure of COX-2 inhibitors, derivatives of dihydropyrimidinone (DHPM) are ideal candidates to be explored as COX-2 inhibitors, due to the ease of synthesis and their versatility to be transformed chemically. In this study, we prepared a novel small library consisting of 288 designed DHPM derivatives by varying the constituent components. The selection criteria of potential candidates for the COX-2 inhibitor of the data bank involve in silico studies via molecular docking investigations, prediction of ADMET and druglikeness, as well as molecular dynamics (MD) simulations. Molecular docking served as the initial step of selection, based on the comparison of grid score, docking pose, and interactions with those of lumiracoxib (LUR) as the original ligand of COX-2. The next criteria of selection were scores obtained from the ADMET and druglikeness by comparing the designed candidates with COX-2 inhibitors that were already marketed. Compound RDUE2 and SDT29 were the most potential candidates, which were further analyzed using the MD simulation. The results of the MD simulation indicated that RDUE2 and SDT29 interacted stably with amino acid residues on the active site of COX-2. The estimation of binding free energy indicated that SDT29 exhibited an inhibitory activity comparable to that of LUR, whereas RDUE2 showed a lower inhibitory activity than that of SDT29 and LUR.


Introduction
Inammation is a response of the immune system to noxious stimuli, such as pathogens, damaged cells, toxic compounds, or irradiation, 1 which is characterized by redness, swelling, heat, pain, and the loss of tissue function. 2This can become chronic if persistent for a long term, which can lead to tissue damage and cell death, causing various kinds of degenerative 3 and neurodegenerative diseases. 4According to WHO, chronic inammation is the most signicant cause of death in the world. 5Diseases associated with chronic inammation are expected to continue to increase over the next 30 years.[8] Generally, anti-inammatory drugs are COX-2 inhibitors, but these drugs cause unwanted side effects, such as the risk of heart and liver diseases. 9,10Therefore, new drug candidates for the treatment of inammation are urgently needed.
The inhibition of COX-2 is the main solution in the treatment of inammation.This enzyme is involved in the conversion of arachidonic acid to prostanoids, which are important substances involved in the occurrence of inammatory processes. 11Therefore, COX-2 can be used as the target for the treatment of inammation.One of the highly potent and selective drugs that inhibits COX-2 is lumiracoxib, which has an IC 50 value of 0.13 mM with an excellent selectivity ratio (IC 50 COX-1/IC 50 COX-2) of 515 in the human whole blood assay. 12his drug has low gastrointestinal side effects; 13 however, recently it has been reported to cause liver damage, leading to its withdrawal in 2007. 14,15Nevertheless, the molecular interaction between this drug and COX-2 serves as a good starting point for designing new inhibitors with different topologies to minimize the possibility of similar side effects.According to Carullo et al., 16 COX-2 inhibitors require two signicant moieties that must be present in their structure, namely, carboxyl group and aromatic ring.These parts can bind with COX-2 to exhibit the desired inhibitory effect.
One of the compounds with diverse activities that can be easily designed according to the required structure is dihydropyrimidinone (DHPM). 17DHPM and their derivatives have been extensively studied and are known to have promising biological activities including anti-inammatory, 18 antioxidant, 19 anticancer, 20 and antimicrobial functions. 21In silico studies proved that DHPM possesses the potential to act as a COX-2 inhibitor compared to Celecoxib. 17Alfayomy et al. 22 revealed that DHPM shows a great IC 50 value compared to Celecoxib.In vivo studies showed that DHPM has low side effects on gastrointestinal ulcers.Based on these ndings, it is assumed that DHPM derivatives have great potential as COX-2 inhibitors.
This research was aimed to discover new COX-2 inhibitors through in silico studies.Candidates with DHPM were rst designed as the main core, and then a docking experiment was performed to understand the inhibitor binding mode to the enzyme.ADMET prediction was used to study the drug-likeness of the candidates.MD simulations and binding free energy calculations were carried out to deeply investigate the interactions that occur between inhibitors and COX-2.This research is intended as a logical approach to the development and nding of potent inhibitors, showing good efficacy and efficiency as COX-2 inhibitors.

Materials
COX-2 crystal complex with lumiracoxib (LUR) (PDB ID: 4OTY) was retrieved from RCSB Protein Data Bank, 23 and 30 already marketed COX-2 inhibitors and 288 prior designed DHPM derivatives were used as candidates.

Small molecular structure preparation
The 3D structure of COX-2 inhibitor candidates was constructed and assigned a protonation state at pH 7.4 using Avogadro. 24he structure was further optimized by a PM7 semi-empirical method implemented using Gaussian16. 25Finally, the addition of charges to the candidates was optimized by the AM1-BCC method using an antechamber program. 26

Molecular docking
A docking experiment was started with the validation of docking parameters that started by ligand and enzyme preparation using the DockPrep features in the Chimera program. 27The addition of charge to the receptor was calculated using the ff14SB method, 28 whereas the non-protein part using the AM1-BCC method. 26The surface was made using the structure of a receptor that does not contain hydrogen atoms using the Write DMS feature in the Chimera program.Furthermore, the surface was used to create spheres using the SPHGEN program. 29The selection of spheres was done using the sphere selector program with a distance of 6.0 Å from the ligand position.Aer the spheres were selected, the simulation box was made using the SHOWBOX program with a radius of 7.0 Å.
The production of the grid was based on the Lennard-Jones model with an attractive and repulsion exponent of 6-9 using the GRID program. 30Finally, the redocking of the original ligand was carried out on the receptor by the exible docking method using the DOCK6 program. 31The pose reproduction was considered successful if the redocking results show the same ligand pose as the original ligand before the docking process and have an RMSD value of #2.0 Å. 32 Then all the ligand candidates were docked to COX-2, and the selection of candidates for further analyses is based on the grid score and poses produced during the docking process.

ADMET analysis
The ADMET properties were predicted using the SwissADME web service 33 and Toxicity Estimation Soware Tool (T.E.S.T) program, 34 and the analysis was performed for 30 marketed COX-2 inhibitors and 288 candidates.Candidate selection criteria were assessed using the selection score (SS), 35 which includes the ADMET parameter and grid score from docking.Candidates that showed a higher SS than the reference (ibuprofen) then proceeded to the MD simulation.

Molecular dynamics simulation
Molecular dynamics (MD) simulation was used to study protein-inhibitor interactions. 36The used ligands and receptors were retrieved from the previously docking process.All simulations were carried out using the AMBER22 program. 37Protein was processed under the ff14SB force eld, 28 and inhibitors were processed under the GAFF force eld. 38The complex in the solution phase was modelled using the TIP3PBOX water model with a radius of 10 Å.A total of four Na + ions were added to neutralize the protein-ligand complex system.
The system minimization was carried out by the steepest descent and conjugate gradient method.The minimization stage consists of three steps using the sander program in AMBER22: (i) minimization of water molecules as a solvent with a force constant of 500 kcal mol −1 Å −1 and hold on the N and C terminal protein (3000 cycles of steepest descent and 2000 cycles of conjugate gradient); (ii) minimization of water molecules and side chain atoms of the complex with a force constant of 500 kcal mol −1 Å −1 and hold on residues on the active side and at the end of the protein (1000 cycles of steepest descent and 2000 cycles of conjugate gradient); and (iii) nally, the minimization of all atoms in the system of 2000 cycles without any restrictions (1000 cycles of steepest descent and 1000 cycles of conjugate gradient).
The heating process was carried out for 200 ps from 0 K to 310 K on the NVT Ensemble.Adjustment of the system density was carried out for 300 ps at 310 K.The equilibration system was carried out with ten consecutive stages using partial restraints of 20.0, 10.0, 5.0, 2.5, 1.0, 0.5, 0.1, 0.05, and 0.001, and without restraints at all; in the residues on the active site, as well as C and N terminals of enzyme for 4500 ps on the NPT ensemble.Then, the MD simulation was carried out for 150 ns

Paper RSC Advances
to nd out the stability of the system during the simulation process.The coordinates are stored every 5 ps.All simulations were carried out using the sander program in AMBER22.

Binding free energy calculation
Binding free energy calculation was used to study the binding affinity of inhibitors on the protein.To calculate the binding free energy, 2000 snapshots obtained from the last 10 ns were used.Molecular Mechanics-Poisson Boltzmann Surface Area (MM-PBSA) and Molecular Mechanics-Generalized Born Surface Area (MM-GBSA) were used in the calculation of binding free energy. 39Thus, the binding free energy (DG bind ) between inhibitors (lig) and COX-2 (rec) was calculated using the following equation: where G com , G rec , and G lig are the free energy from complexes, receptors, and ligands, respectively.DG bind gives the gas phase enthalpy (DE MM ), solvation-free energy (DG sol ), and conformational entropy (−TDS) in the binding of ligands.DE MM is the amount of internal energy/DE int (bond, angle, and dihedral energies), electrostatic interactions/DE ele , and van der Waals interactions/DE vdw .DG sol is the number of polar contributions (DG GB/PB ) and nonpolar contributions (DG SA ).
The polar contribution was calculated by the PB and GB approach using programs available in AMBER22. 37We used 80 for the exterior dielectric constant and 1 for the solute dielectric constant.The nonpolar part of desolvation was estimated from the solvent accessible surface area (SASA) using the LCPO method with a water probe radius of 1.4 Å: DG SA = 0.0072 × DSASA.Here, changes in conformational entropy −TDS are not considered due to expensive computational costs and low prediction accuracy. 40

Trajectory analysis
Detailed information on protein-ligand binding was obtained by studying the contribution of each amino acid residue to the interaction energy between the inhibitor and COX-2 using the MM-PBSA and MM-GBSA methods in AMBER22. 41The binding contribution of each residue includes four things: van der Waals (DE vdw ), electrostatic (DE ele ), polar solvation (DG GB/PB ), and nonpolar solvation (DG SA ), without considering the entropy contribution.
In addition, trajectory analysis also includes RMSF, RoG, and SASA analyses to study the stability of each system. 42Trajectory analysis was performed at 2000 frames from the last 10 ns of the MD trajectory.

Design of inhibitors
Generally, COX-2 inhibitors possess a carboxyl group and an aromatic ring in their structure. 11,43The carboxyl group of COX-2 inhibitors interacts with the enzyme through hydrogen bonds, while the aromatic ring interacts at the non-polar part of the COX-2 active site.Such structure is found in many anti-inammatory drugs, such as lumiracoxib, diclofenac, etodolac, urbiprofen, and ketoprofen. 44HPM is a class of compounds that are easy to be synthesized and chemically transformed to meet the structural requirements needed for COX-2 inhibitors.The Biginelli reaction is usually used in the synthesis of DHPM derivatives by mixing the reaction components, namely, derivatives of aldehyde (A), 1,3-dicarbonyl (B), and urea (C) in a reaction vessel. 45or the aldehyde component, several substituted benzaldehydes, heteroaromatics, and cinnamaldehyde derivatives were used.In the carbonyl component, we used 3 aliphatic 1,3diketones.For the urea component, simple urea and cyclic urea derivatives were employed.A small library was created by enumerating these three components, resulting in a total of 288 compounds, including their stereoisomers.The rst letter of the compound name refers to absolute conguration information (ESI, Table S1 †).The core structure of the candidates is displayed in Table 1.

Docking analysis
Docking was used to study the interaction mode of the inhibitor on COX-2.The docking of the original ligand, lumiracoxib (LUR), gave a grid score of −64.16 kcal mol −1 .The carboxyl group of LUR built two hydrogen bonds with Tyr354 and Ser499.The 3D and 2D visualization of the interaction between LUR and COX-2 is presented in Fig. 1.

Component variations
Main core of candidates Then, we docked 288 designed candidates and 30 marketed COX-2 inhibitors.The information obtained from the results comprised grid score, docking pose, and interaction of candidates with COX-2.To maximize prediction accuracy, the docking poses of candidates should overlap with the cognate ligand (LUR) and the candidates with similar interaction modes to LUR are expected to show better or at least equivalent activity to LUR. 46,47 The results showed that only one candidate had a lower grid score than that of LUR.However, its docking pose did not have same interaction and did not overlap with LUR, so it was not selected for further analyses.
The docking results of LUR and 14 selected candidates are shown in Table 2. Generally, the carboxyl group of the candidates forms hydrogen bonds with residues Tyr354 and Ser499.Based on these data, it is known that DHPM synthesized from the aldehyde component with a phenoxyacetate moiety predominantly occupies the top rank.The selection process was then followed by the ADMET and druglikeness prediction.

Selection of candidates
The selection of candidates was based on ADMET and docking parameters.Based on these parameters, Castro-González et al. 35 created an assessment model called the selection score (SS) to sort drug candidates by comparing the parameter scores of the candidate with the parameter score of ibuprofen as drug references.Ibuprofen was used as a reference because it is well known for its activity as an COX-2 inhibitor and in general has been used by the public as an anti-inammatory drug. 48The results of the selection score analysis of 14 selected candidates are displayed in Fig. 2.

Molecular dynamics simulation
Structural exibility and stability.In general, molecular dynamics simulation was carried out on a protein-ligand system to explore its stability and energy.The MD simulation was run for 150 ns on both the apo-COX-2 system and the protein-ligand system (COX-2:LUR, COX-2:RDUE2, and COX-2:SDT29).The complex stability was assessed based on the root mean square deviation (RMSD) of the atoms in the backbone and ligands.
Fig. 4 describes the RMSD stability of the backbone and ligand during simulation, and it was achieved during the last 10 ns of the simulation, ranging from 0.17 to 0.24 nm and 0.02 to 0.10 nm, respectively.The apo-COX-2 system showed the largest RMSD backbone, which indicated that the free COX-2 structure tends to have a greater structural change than that in a complex system.
However, the parameters of stability were determined not only by RMSD value but also by free energy.Therefore, the next parameter of complex stability to be determined was energy stability during simulation.This parameter was obtained from the MM-GBSA and MM-PBSA approach, including the energy of non-polar (van der Waals) and polar interactions (electrostatic interaction and hydrogen bonding), as shown in Fig. 5.
The total (nonpolar and polar) energy analyzed using MM-GBSA and MM-PBSA approaches was stable, which declared that all systems were relatively stable, and they did not require an additional longer time of simulation.In addition, all three systems achieved their stability at about the same time.
The achievement of stability as shown in the analysis of RMSD and total energy was then continued with trajectory analysis, which was carried out for the last 10 ns of the simulation (140-150 ns).It comprises root mean square uctuation (RMSF), radius of gyration (RoG), and solvent accessible surface area (SASA).
RMSF analysis was used to determine the system exibility through uctuations of amino acid residues. 49Fig. 6A shows that all systems uctuate at residues of 1-2, 18-22, 41-53, 66-67, and 552.Specically, residue numbers 183-184 uctuate in apo-COX-2 and COX-2: SDT29.Several uctuations only occur in COX-2: LUR, namely, residues 227, 343, and 511, while uctuations in residue numbers 339 and 523 only occur in COX-2:SDT29.Overall, the uctuations that occurred in each system are residues that are quite far from the active site of COX-2.This observation gave us information that the active site of the enzyme has low exibility, which allows for a more stable interaction between ligands and COX-2.
The compactness of systems was assessed by the RoG of each system, as shown in Fig. 6B.It shows that each system has relatively the same RoG but with a quite small deviation, so that it can be assumed that the entire system is stable and compact.Fig. 6B shows that the RoG of the system from smallest to largest is: COX-2:LUR, COX-2:SDT29, COX-2: RDUE2, and apo-COX-2.According to Khan et al., 50 the lower the RoG of a system, the more compact the system is.The system with a ligand on the active site of COX-2 possessed a lower RoG than the apo system.This proves that ligands that interact with COX-2 can affect the structure and compactness of the overall structure of the enzyme.
Solvent accessible surface area (SASA) is one of the determining factors in studying the folding and stability of enzyme structures. 51Enzyme structures that experience an unfolded state tend to have a higher SASA than those with stably folded structures. 52In this study, SASA analysis was performed on the active site and the entire enzyme.The results indicated that for both the active site and the whole enzyme, apo-COX-2 had the highest SASA, while COX-2:LUR provided the smallest SASA.This indicates that the presence of a ligand on the active site of COX-2 results in a more stable complex structure than the structure of the enzyme without any ligand.The presence of candidates certainly also affects the stability of the folding of COX-2, which can be seen by its SASA.The average SASA of each system is shown in Fig. 6C and D.
Binding free energy calculation.The stability and affinity of each complex are represented by the binding free energy (DG bind ). 53MM-PBSA and MM-GBSA were used to estimate DG bind when ligands interact with COX-2. 39The last 10 ns trajectory of the simulation was used to calculate DG bind , under consideration that the stability of each system had been achieved.It means that there are no signicant uctuations in this trajectory, so that it is useful for increasing the efficiency of calculation when conducting analysis. 42DG bind for each system is shown in Table 3.
The results revealed that COX-2:RDUE has the largest DG bind value compared to others.While COX-2:SDT29 exhibits a similar DG bind value to COX-2:LUR used as a reference.This indicated that SDT29 has an activity equivalent to LUR and better than RDUE2.
Analysis of energy decomposition.Energy decomposition analysis was used to nd out which contributions are given by residues for the stability of the ligand-enzyme complex.Based on the results, only a small portion of residues greatly contribute to binding free energy.In addition, there are no residues that contribute to unfavorable interactions such as repulsion between ligands and enzymes.The graph of energy decomposition analysis is shown in Fig. 7.
Energy decomposition analysis shows that Tyr354 and Ser499 are residues that gave the greatest contribution to the binding free energy of the three complexes.However, in COX-2:SDT29, Val492 also gave a signicant contribution to the stability of the complex.The presence of this additional contribution causes COX-2:SDT29 to have good stability and a lower binding free energy compared to COX-2:RDUE2.It is also conrmed that SDT29 inhibition activity towards COX-2 is better and comparable to LUR than RDUE2.The analysis of the type of contribution for each residue based on decomposition energy is shown in Fig. 8A and Tables S4, S5. † It is shown that Tyr354 and Ser499 are two benecial residues for the stability of the three complexes, by contributing to electrostatic energy.This is in line with the results of molecular docking analysis, where ligands form hydrogen bonds with COX-2 through these two residues.However, the energy of the polar solvation of these residues is unfavorable for the stability of the complex.In the complex of COX-2:LUR and COX-2:RDUE2, Val492 provided electrostatic energy that is not benecial for complex stability.However, Val492 shows bene-cial electrostatic energy in COX-2:SDT29.
In addition to the electrostatic energy, several other residues make a quite good contribution to the stability of the complex through the van der Waals interaction.The residues include Tyr317, Val318, Leu321, Ser322, Phe487, Val492, Gly495 and Ala496.These residues interact with non-polar parts of ligands, such as aromatic rings in LUR, RDUE2, and SDT29.
Residue interaction.The interaction of ligands with residues around the active side of COX-2 was studied with contact analysis performed on the last 10 ns trajectory.The contact analysis process was focused on amino acid residues with a radius of 3.5 Å from the active side of the enzyme.In general, there are two types of interactions between ligands and residues on each complex, namely, hydrogen bonds and van der Waals interactions.Contact analysis is focused on the percentage of interactions possessing a fraction >1%.The fraction can be calculated using eqn (5), which shows the percentage of interaction (fraction) of the number of frames for each interaction (N fra ) divided by the total number of frames (N tot ) during the simulation of 2000 frames. 54action ¼ Hydrogen bond interaction between residues and ligands showed high fraction and was assumed to contribute most of the complex stability in each system.Fig. 8A shows that during the last 10 ns of the simulation time, on average there are two hydrogen bonds that occur in the COX-2:LUR and COX-2:RDUE2 systems.However, the occurrence of hydrogen bonds in the COX-2:LUR system is considered more consistent than that in the COX-2:RDUE2.However, on average, there are three bonds that occur in the COX-2:SDT29 system.This shows that the stability of COX-2:SDT29 is better than that of COX-2:RDUE2.This data showed that Ser499 and Tyr354 are the two main residues playing a major role in hydrogen bonds in each complex.Besides both interactions, in complex COX-2:SDT29, an additional hydrogen bond interaction with Val492 in lower fractions was observed.
In addition to the hydrogen bonds, complex stability was supported by van der Waals interactions between ligands and COX-2, as shown in Fig. 8A.Contact analysis informed that LUR, RDUE2, and SDT29 also interact with COX-2 via van der Waals interactions through Ser499 in high fractions.Besides Ser499,   RDUE2 shows van der Waals interactions in high fractions with Leu321, while SDT29 shows van der Waals interactions with Val318 and Tyr324.In addition to the interaction of van der Waals in high fractions, there are quite several other van der Waals interactions in lower fractions, which were expected to support the stability of the complex.Briey, all these observed interactions between ligands and receptors supported the inhibitory activity against COX-2.

Conclusions
A series of DHPM derivatives have been designed as COX-2 inhibitors.Their inhibitory activities were deduced from the results of molecular docking, ADMET prediction, drug-likeness, and molecular dynamics simulation.The selection of the potential candidates was based on the criteria reported by Castro-González et al. 35 The structure of COX-2 in the complex with lumiracoxib was retrieved from the Protein Data Bank (PDB ID: 4OTY), while 288 candidates were designed of dihydropyrimidinone scaffolds as products of the Biginelli reaction, and 30 drugs known as COX-2 inhibitors were used as reference in drug-likeness analysis.Thirty two of 288 candidates showed better grid score and similar binding mode to LUR, which were then taken for the selection of drug-likeness, which comprised pharmacokinetics and toxicity.The selection scores obtained were then used to determine the potential candidates.Based on the in silico study, two candidates, namely, RDUE2 and SDT29 are predicted to possess potential inhibitory activities, which are comparable to LUR.Amino acid residues playing a pivotal role in the inhibition activity are Tyr354 and Ser499 via the formation of hydrogen bonding with the carboxyl group of the ligand.
For molecular dynamics simulations, apo-COX-2, COX-2:LUR, COX-2:RDUE2, and COX-2:SDT29 systems were used to study the stability and energy of each system.The estimation of binding free energy, energy decomposition analysis, and interaction analysis between residues and ligands indicate that RDUE2 and SDT29 have promising potential as COX-2 inhibitors compared to LUR.The binding free energy of SDT29 was smaller, and the residue number of Ser499, Tyr354, and Val492 contributed to the stability of COX-2:SDT29.However, the binding free energy of COX-2:RDUE2 was greater, and only Ser499 and Tyr354 contribute to their stability.The energy decomposition analysis showed that the main residue contributes to complex stability via electrostatic interaction.The results indicated that the hydrogen bonds between candidates and COX-2 are very useful for increasing the inhibitory activity.These results are following our assumptions when designing new COX-2 inhibitors based on DHPM, and the carboxyl group is an important part because it endows the ligands with the ability to interact with COX-2 through hydrogen bonds, causing the expected inhibitory activities.

Fig. 2
Fig.2Graph of SS analysis of the selected candidates and ibuprofen.

Fig. 5
Fig. 5 Stability of (A) total nonpolar energy (DE vdw + DG SA ), (B) total polar energy (DE ele + DG GB ) based on the MM-GBSA approach, (C) total nonpolar energy (DE vdw + DG SA ), and (D) total polar energy (DE ele + DG PB ) based on the MM-PBSA approach.

Fig. 6
Fig. 6 Trajectory analysis of each system: (A) RMSF, (B) RoG, (C) SASA on the active site, and (D) SASA on the entire enzyme.

Fig. 7
Fig.7Energy decomposition analysis of residues which greatly contribute to the stability of the complex.

Fig. 8
Fig. 8 Contribution for each residue based on the decomposition energy for COX-2 (A) and residue interactions with (B) LUR, (C) RDUE2, and (D) SDT29.

Table 2
Docking results of LUR and 14 selected candidates